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Overview 


The objective of this research is to study aeroacoustic noise generation in a supersonic 
round jet. In particular we want to understand the effect of turbulence structure on the 
noise without numerically compromising the turbulence itself. This means that direct 
numerical simulations are needed. In order to use DNS at high enough Reynolds numbers 
to get sufficient turbulence structure we have decided to solve the temporal jet problem, 
using periodicity in the direction of the jet axis. Physically this means that turbulent 
structures in the jet are repeated in successive downstream cells instead of being gradually 
modified downstream into a jet plume. Therefore in order to answer some questions about 
the turbulence we will partially compromise the overall structure of the jet. 

This report will be presented in two chapters. The first chapter is divided into two 
sections. The first section describes some work on the linear stability of a supersonic round 
jet and the implications of this for the jet noise problem. In the second section we will 
present preliminary work done using a TVD numerical scheme on the CM5 at the Univer- 
sity of Minnesota. This work is only two.-dimensional (plane) but shows very interesting 
results, including weak shock waves. However this is a non-viscous computation and the 
method resolves the shocks by adding extra numerical dissipation where the gradients are 
large. One wonders whether the extra dissipation would influence small turbulent struc- 
tures like small intense vortices. Nevertheless we are still pursuing a three-dimensional 
version using this method for the purpose of comparing with other methods, and perhaps 
answering some of our doubts. 

The second chapter is an extensive discussion of preliminary numerical work using the 
spectral method which we prefer to use to solve the compressible Navier-Stokes equations 
to study turbulent jet flows. The method uses Fourier expansions in the azimuthal and 
streamwise direction and a 1-D B-spline basis representation in the radial direction. The 
B-spline basis is locally supported and this ensures block diagonal matrix equations which 
are solved in O(N) steps. A very accurate highly resolved direct numerical simulation 
(DNS) of a turbulent jet flow is expected. This is a modification of a boundary layer code 
developed by Bob Moser. 
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Chapter 1 : Preliminary Studies of Inviscid Compres sible Jets 
Section 1 : The stability of a supersonic round jet 

Our purpose in doing this stability problem is to provide initial conditions for the 
nonlinear computations. We want to start with the most unstable disturbance since this is 
what would be found in the natural problem. This problem has been studied by Tam and 
Hu(1989) however they didn’t present the information which we require. As a result of our 
computations we have found some new results which have not been previously reported. 

The basic stability equation was obtained by linearizing the equations for inviscid 
compressible flow in cylindrical coordinates, taking small perturbations from an axially 
symmetric basic state. One can reduce this to a single second order equation for the 
pressure perturbation, which is expressed in separation of variables form as 

p = p(r)+p{r)e^ t - ikz - ime] 

where r is the radial coordinate. The equation for p is 

d 2 p 1 1 2 kdu/dr 1 dp\dp / (oi-ufc) 2 fc2 _ - Q 

dr^^Kr^ uj - uk p dr ) dr V a 2 7-2 2 

where u{r) is the prescribed axial velocity profile, p{r) and p(r) the prescribed pressure 
and density of the base state, and a 2 (= jRf) is the square of the sound speed. This 
equation agrees with Tam and Hu. This equation is to be solved with boundary conditions 
which require the solution to be bounded at r = 0 and to possess only outgoing waves at 
infinity. It is an eigenvalue problem for complex ui when k and m are given. The flow is 

unstable if Imu; < 0. 

The velocity profile was taken in the same form as in Tam and Hu as a “half-Gaussian” 
function which is given by 

■u = uj r < h 

u = Uj exp ^ — ln2( — — — ) ) r > h 
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Half - Gaussian jet profile 



which is sketched here in order to define the parameters b and h. We define the jet radius 
to be rj = b + h and all computations were done with b/rj = .1, h/rj = .9. We have taken 
the temperature profile to be of a form similar to the velocity 


T = Tj r < h 

T = r oo + (T i -r oo )exp(-ln2(^) 2 ) r>h 
and have taken p(r) = constant across the jet. Then p is related to T by 


Poo = pRT 


and 

® = a^T/Too- 

Sometimes in problems like this the temperature is related to the velocity profile by the 
Crocco-Buseman formula. But this is a viscous dominated profile which results when the 
Prandtl number is unity, and it doesnt make a lot of sense for high Reynolds number flows 
such as this. In all the work done by us so far we have taken Tj = so the temperature 
is constant across the jet. 

We have solved this eigenvalue problem by a method described in a book by Betchov 
and Criminale. We integrate inward from a large radius, using an asymptotic formula 
to identify outgoing waves, and we integrate outward from the origin. Where the two 
computations meet we must have p and dp/dr continuous. This does not occur unless 
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u! has the proper value. We iterate on u> using a Newton-Raphson procedure until our 
continuity requirement is attained. The biggest problem was that the procedure doesn’t 
converge unless one starts near the proper value. This means that is was necessary to 
change parameters by small increments, starting from a known result. 

Results have been obtained for jet Mach numbers of Mj = 2 and 2.5. The result for 
Mj = 2 is shown in a figure at the end of this section. Here we present curves of growth 
rate versus k for a number of values of m (the azimuthal wavenumber). What is observed is 
that the maximum growth rate for given m increases as m increases, reaching a maximum 
at m = 4 and then it decreases again as m becomes greater that 4. This means that of all 
azimuthal wavenumbers the most unstable wave is m = 4 (and also m = — 4 since only m 2 
occurs in the equation). For Mj = 2.5 the maximum is at m = ±5. 

A wave exp(— ikz — imd ) is interpreted as a helical wave on the surface of the jet, 
with the wave making an angle tan -1 (ra/fc) with the axis of the jet (m = 0 is an axial 
wave). Therefore there are two helical waves of definite angle, one left-handed and one 
right-handed, which are the most unstable. This is similar to results of Sandham and 
Reynolds (1990) on the stability of a compressible mixing layer. They find two oblique 
waves with maximum growth rate which make equal but opposite angles with the flow 
direction. 

There are some implications of this for our proposed nonlinear computation. Sandham 
and Reynolds (1991) in a second paper find a single oblique unstable plane wave rolls up 
into an oblique vortex. We expect that a single unstable helical wave will roll up into a 
helical vortex. Since our helical waves propagate downstream with a speed which is about 
.66 times the jet speed, we expect a similar speed for the helical vortex. A propagating 
helix will appear like a rotating helix-as a barber pole does. This is an interesting result 
because there are observed jet modes, called spinning modes, where the noise field rotates 
around the jet axis, so that to an observer on one side the sound appears to have a time 
periodicity (Wesley and Woolley, 1975). 

Sandham and Reynolds found that the angle of the most unstable oblique wave de- 
pends on the convective Mach number in such a way that M c cos a ~ .6, so that the 
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component of the Mach number perpendicular to the wave stays constant (and subsonic). 
We find approximately the same result for the helical waves, based on only two Mach 
numbers. In the Sandham-Reynolds nonlinear computations they find that because the 
normal Mach number is subsonic the oblique vortex forms without shock waves (unlike 
their strictly two-dimensional computations which do have shocks). This is possibly very 
significant for our proposed computations. It means that helical vortices could form with- 
out shocks. This is important when we use the B-spline/spectral method because any 
shocks will have to be resolved by viscosity alone which could restrict our Reynolds num- 
ber range. 
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Section 2 : A Study of 2-D, Unsteady, Inviscid, Compressible Jets Using 

a Second Order TVD Scheme 

Computations of a two-dimensional supersonic temporal jet were made as a prelim- 
inary study of the non-reflecting boundary conditions and to gain some experience with 
the temporal jet configuration. This was done using an efficient, parallel TVD 5-7 scheme 
in a data-parallel environment to run on massively parallel processing (MPP) computers. 
All the calculations are carried out on Thinking Machines’ CM-5. It is hoped that these 
calculations will help discern the intricacies of 2-D compressible jet flows. An attempt is 
made to resolve the turbulent scales(vortex structures) using very fine meshes. No turbu- 
lence model is used. Since these are inviscid calculations the numerical dissipation of the 
scheme is used to approximate viscous dissipation effects of real fluids. In this study a 
second order explicit non-MUSCL upwind TVD algorithm as proposed by Yee 7 is used. 

Initially a 2-D version of the half-gaussian jet profile described in the previous section 
is used. Anti-symmetric disturbances were imposed on the initial jet profile since they 
were found to be more unstable. This ensures the break down of the shear layer into small 
vortical structures. 

The computational domain was discretized into a very fine mesh of 512*512 grid cells. 
The physical dimensions of the domain are twice the wavelength of the most unstable wave 
(found from linear stability analysis) in the streamwise direction and 3 to 4 jet diameters in 
the transverse direction. Grid stretching is used for a layer of cells (one jet diameter wide) 
near the boundary in order to weaken any outgoing waves. At the transverse boundaries 
the approximate non-reflecting boundary condition 10 p' ± pav 1 = 0 was used. The primed 
quantities are perturbations from mean flow, and a is the speed of sound. 

Simulations for a fully expanded cold jet at centerline Mach number of 2.0 are pre- 
sented. 

Highly resolved flow solutions capturing a wide range of turbulent scales are achieved. 
Figures 1-3 show the various stages of temporal jet development. The Mach contour 
and pressure contours show distinct development of shocks which were absent initially. 
The isobars also show alternating high and low pressure regions corresponding to vortical 
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structures. From the vorticity contours we notice that the initial structure of two vortex 
sheets breaks down into smaller vortices. 

It is evident from the plots that the non-reflecting conditions work fairly well. 
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Chapter 2 : Direct Numerical Simulation of Unsteady, Compressible Round Jets 


Introduction 


Jet noise is a major concern in the design of a supersonic transport (SST) aircraft. 
Studies by various researchers 2-4 lead to aerodynamic noise as a major contributor to 
jet noise. Some of these studies indicate that most of the aerodynamic jet noise due to 
turbulent mixing occurs when there is a rapid variation in turbulent structure, i.e. rapidly 
growing or decaying vortices. 

The object of this work is to obtain highly accurate flow solutions of a turbulent 
round jet. These solutions are expected to help understand the various turbulent scales 
and mechanisms of turbulence generation in the evolution of a compressible round jet. We 
hope to use these accurate flow solutions to estimate acoustic radiation in the near-field 
region. Also the data generated can be used to compute various turbulence quantities such 
as mean velocities, turbulent stresses, etc. which may aid in turbulence modeling. 

We simulate a compressible round jet by using Fourier expansions in the azimuthal 
and streamwise direction and a 1-D B-spline basis representation in the radial direction. 
This is an efficient and accurate way to separate out the 0 and z variables, leaving partial 
differential equations(PDEs) depending on r and t only. By using a 1-D B-spline basis 
and a Galerkin approximation we can reduce this set of PDEs to ordinary differential 
equations(ODEs) in time. This is solved by a 3 rd order Runge-Kutta time marching- 
scheme. The present study uses a spectral method developed by Moser et a l. 7 

We consider the temporal jet problem for two reasons: one since this configuration 
allows for the application of highly accurate spectral methods and the other because the 
dynamics of the temporal jet is not greatly different from that of a spatially evolving jet. 
The spectral accuracy helps capture smaller turbulent scales. 

Governing Equations 

The compressible Navier-Stokes equations written in cylindrical coordinates in non- 
dimensional form are, 
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where, 


_d_ _ d_ JL-1JL d _ d 

dx 1 dr ’ c?X 2 r dO ’ 9x3 9^: 


a = p> m k = P u k i -Re = Pr = 

i?e, Pr are the Reynolds number and Prandtl number, is the specific heat at 
constant pressure, and $ is the viscous dissipation (see Appendix). 


B-Spline Representation and Galerkin Formulation 


The flow variables are expanded using Fourier sums in the two periodic directions, viz. 
the azimuthal (6) and the axial (z) directions. In the non-periodic or radial direction (r) 
we use 1-Dimensional B-splines as interpolating functions. 1 B-splines have local support 
and hence lead to sparse block diagonal matrices which can be efficiently stored and solved. 
Application of boundary conditions is similar in ease to a finite element method (FEM). 

B-splines of order n are piecewise polynomials of degree n having n — 1 continuous 
derivatives. Since they have a high degree of continuity derivatives of quantities (like 
velocity to get vorticity) can be smoothly and accurately represented. 

They have 1 degree of freedom(d.o.f) per interval unlike finite element(FE) basis which 
can have as many d.o.f as the order of the polynomial. We use B-spline bases because 
higher order FE bases may resolve waves of wavelength smaller than twice the interval 
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length due to the increased d.o.f. But this is beyond the Nyquist cut-off for the interval 
size and so these increased d.o.f do not improve the resolution (i.e. the smallest accurately 
represented scale) of the solution only thing that gets better is the accuracy (convergence 
rate is increased). Higher order B-splines on the other hand not only have better accuracy 
but also have better resolution of scales per d.o.f. 

The Galerkin method using B-splines as basis functions has been used previously by 
some researchers. 5-7 By using a Galerkin approximation we can approximate the set of 
PDEs as a set of ODEs in time. 

One can use a B-spline basis to represent the desired function f(x) as, 

OO 

/(*) = ^ a 3 b ]( X ) ( 4 ) 

j=z~00 

where, 

b’j(x) is a n th order B-spline coefficient, a,j is the value of function at knot j. 

Using a Galerkin approximation we can write, ( b k is the weight function) 

j b k f(x)dx = 'Y^a j J bjb k dx (5) 

j 

Also, the derivative of f(x) can be written as, 


/'(*) = £> (.'(*) (6) 

j 

But here the order of the polynomial has reduced by 1 so in order to keep the degree 
of polynomial the same we approximate the derivative as, 


f{x)*ig(x) =► J2 a i b j( x ) ~ 


Again the Galerkin approximation is written as, 


J b k f'(x)dx ss ^ a,j J bj b k dx 


(7) 


( 8 ) 
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Any non-linear terms are handled in a similar way, i.e. if 

h = f g 

h - 22 dj bj(x) R 2 22 a i Ck bk 

3 

'y > dj J bj bi dx ~ y ^ dj Ck / bj bk b[ dx (9) 

j j,k 

The matrices (terms with the integral) on the right hand side of equations (5) and 
(8) are called the mass and derivative matrices respectively. All different combinations of 
such derivatives and other terms are computed as matrices too. These are calculated only 
once using a Gaussian quadrature and stored as opposed to a regular finite element (FE) 
basis where they can be calculated on the fly when required. 

Writing the governing equations in the Galerkin form using the B-spline representa- 
tion yields a number of matrices similar to the mass and derivative matrices, non-linear 
advection terms yield matrices similar to the non-linear matrix obtained in equation (9). 

In Galerkin form the continuity equation can be written for any kg and k z as (kg and 
k z are wave numbers in the 8 and z directions respectively), 


v - '' f , , , ( f h bj (r bk)' bi 

2 / bi bi rdr = 2_^ &i <J j( m r k / — — rdr 

i Jr „• i. i Jr 


j , k , l 

d 

+ 


bi bj bk b\ 


rdr -I- j-m Zk J bi bj b k bi rdr^j 


where, 

a = ]T>(z,M) b k (r), m e = y^mg k ( z ,9,t) b k (r) , 

k k 

m r = 22 mr k (*, M r ) ’ = bk(r) 

k k 


The Fourier terms (e lke e tk *) are included in the coefficients of the variables. So 
cr k (z,6,t) = a k (t) Y, kg Hk z e lkf>e+ik > z and so on. 

Since we are using 1-D B-splines all the integrals are line integrals over the radius. 
These are computed exactly using Gaussian quadratures (doing integrals exactly takes care 
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of aliasing). The derivatives in the 9 and 2 are taken by Fast Fourier Transforming (FFT) 
into wave space and multiplying by the appropriate wave numbers kg and k z , and then an 
inverse FFT is applied to bring it back to physical space. 

The momentum and energy equations can be written in a similar manner. 



Fiq 1 Quadratic B-Splines 


Non-Reflecting Boundary 



Fig 2 Computational Domain 

Numerical Formulation 

Writing the flow equations in a manner as discussed above results in a linear system 
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of coupled equations to solve simultaneously at each time step. Since these B-splines (of 
order n) have local support on n + 1 knot(node) intervals (see fig.l) we get a 2n + 1 block 
banded matrix system. 


M f = R (11) 

where, 

M is the resulting mass matrix, / is the column vector of nodal values we solve for, 
and R is a column matrix resulting from the RHS of the governing equations. 

Time integration is carried out by a 3 rd order Runge-Kutta scheme. 

Regularity Requirements 

In the cylindrical coordinate system the origin (r — 0) is source of concern since some 
of the functions do not remain analytic as they have a V in the denominator. From a 
mathematical point of view the flow variables should be single valued and finite. To enforce 
this the polynomial expansion functions must satisfy some regularity requirements. 8 

The z-component of the velocity should be represented as, 


u z (r; m, k ) = a(m, k ) P z (r 2 ; m, k ) e l me e l kz , m = all integers , 

Pz{ 0 ; m, k) = 1 , 


( 12 ) 


where P z (r 2 ; m, k ) is a polynomial in r 2 . 

Scalars and z-components of all vectors should be represented in a similar manner. The 
9 and r-components of the vectors are dependent on each other and should be represented 


as, 


u r (r; m, k ) 
u e (r- m,k ) 
c(m, k) 
c(m, k) 
P r ( 0; m, k) 


= b(m, k) r 1 ™ 1 - 1 P r (r 2 ; m, k) e l me e i kz 
= c(m,k) r^ -1 Pg(r 2 ; m, k) e l m e e l kz 
= ib(m,k) for m > 1, 

= —ib(m,k) form < — 1, 

= P e { 0; m, k ) = 1, 


(13) 
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for to = 0 b(m, k ) and c(m, fc) are unrelated. 

Enforcing these conditions gives rise to a set of constraint equations which then replace 
some rows in the mass matrix and suitably modify the RHS vector R. 

We can obtain the constraint equations as follows, for a quadratic B-Spline and any 
scalar u , 


u(r; m,k) = a(m,k ) r ^ P(r 2 \ m,k ) e l m 9 e l kz , m = all integers , 

Only non-zero derivatives allowed are |m| +2j Vj > 0 

For quadratic B-splines (see hg.l), at r = 0 , 3 splines have support at the origin. 
Order of non-zero derivatives for the splines is, 

Spline 1 =>• 0, 1, 2 
Spline 2 => 1, 2 
Spline 3 => 2 

All other splines are zero at r = 0. 

Consider spline expansions as Y2i a i > where are the spline coefficients. We have 
in any interval, 


^ ] bi — 1 bi + &2 + 63 — 1 
i 

b [ + 6' 2 = 0 =» 6^ = — 62 

Now consider different values of the azimuthal wave number m, 
for m = 0 

Only non-zero derivatives allowed are 0, 2, . . . even powers, so all odd powered deriva- 
tives should be forced to zero. At r = 0 splines having non-zero first derivative are splines 
1 and 2, so, 


a 1 fq + 0262 — 0 


(14) 
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for m = 1 


Only non-zero derivatives allowed are 1,3,... odd powers, so all even powered deriva- 
tives should be forced to zero. At r = 0 splines having non-zero 0 th and 2 nd derivatives 
are splines 1 and splines 1, 2, and 3 respectively, giving constraints, 


a i = 0, and 
a 2^2 ~b a 3^3 = 0 


(15) 


Similarly, 
for rn — 2 

Only non-zero derivatives allowed are 2, 4, . . . even powers, giving 

Oi = 0, 02 = 0 (16) 


and 

for rn = 3 

Only non-zero derivatives allowed are 3, 5, ... odd powers, so 

«i = a 2 = a 3 = 0 (17) 

or, all splines should be zero. For all other values of m we get the same result as (17). 
Now equations (14-17) are the constraint equations which will replace the appropriate 
rows in the matrix M and vector R. 

Before replacing rows which actually contain the physics of the flow care has to be 
taken to see that no information is lost. So, we need to take a linear combination of those 
rows to be replaced and add them to the remaining rows in a particular manner. For this 
the null space of the constraint equations has to be computed and the eigenvectors of this 
space are to be used for the linear combination. This is done as shown below. 

for o]2gl = 7 and m = 1 

The first set of rows in the mass matrix are written as, 
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The rows with *’s are the unconstrained rows and the b ’• are the bspline derivatives. 
Now to choose a set of null vectors Xj such that, 

m'Xj = 0 


(18) 


Let us choose 




0 

0 

1 

0 


0 

0 

1 

£12 

£21 

0 

0 

V^13 

X22 


0 

0 

0 

0 

0 

1 

3-31 


\ 


So we can compute the null vectors using equation 18. 

Now the mass matrix M is modified by using the null vectors for linear combinations 


as, 


row 2 = row 2 (original) + in * row 2 + * row 5 + x 13 * row 7 

row± = rou)4(original) + x 2 i * row 5 + x 22 * row 7 

rowe = rowQ(original) + x 2i * row 7 
This has to be done for all values of the azimuthal wave number m. 

This procedure can similarly be applied to all the flow vectors and scalars appearing 
in the vector /. The derivation of constraints is similar for higher order splines. 

CFL and Modal Reduction 

Another concern arising due to a cylindrical mesh is that near the origin the aspect 
ratio of the cells gets very high, so we have to reduce the number of azimuthal modes 
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near the origin to maintain a good CFL number. This we will call mode suppression. So 
we effectively reduce the number of azimuthal modes to 1 near the origin and increase it 
successively such that at the outer boundary we have all the azimuthal modes. 

This reduces the accuracy but helps us increase the time-step, dt, by a significant 
amount. Another advantage of doing this is that the computational domain in Fourier 
space is significantly smaller thus allowing faster computations, i.e. For lower azimuthal 
wave number x, the computations have to carried out for more radial points, but as x 
increases fewer radial points have to be considered. Another way to look at it would be 
near the origin fewer x loops need to be considered. 

This is to be implemented by making the upper bound of the rc-loops to be dependent 
on y, i.e. x= 1 to nx-modified[y] if the y-loop is the outer loop, or making the lower bound 
of the y-loop dependent on x if the x-loop is the outer loop, i.e. y = ny-min[x] to ny. 

Optimizations 


The code was optimized thus, 

• Moved all major computations (rhs of governing equations) to separate subroutines 
thus reducing the load on the main program and enabling it to be unrolled and vec- 
torized more efficiently. 

• Manually unrolled the quadrature sum loops. This helped in decreasing the number 
of operations by eliminating a lot of repetitive computations due to proper regroup- 
ing. Similar and symmetric matrix multipliers were regrouped together thus reducing 
redundant calculations. 

• Vectorizing over larger loops. 

Following the above procedure has helped reduce the time taken per mode per Runge- 
Kutta step considerably. 


Boundary conditions 

Ideally, the computational domain should mimic the physical domain by including 
all free space and only having physical boundaries. But this is not possible and the 


21 



computational domain has to be finite. So we need artificial boundaries to limit the 
computational domain. But we want these artificial boundaries to be invisible to the 
flow field so that vortices and other waves can pass through these boundaries and leave 
the domain without giving rise to spurious reflection waves. So we need non-reflecting 
boundary conditions (NRBCs). 

In solving a temporally evolving jet the inflow-outflow boundaries are made periodic 
(see fig 2). This might not be very accurate but serves the purpose of studying turbulence. 
Also it takes care of the inflow-outflow boundary conditions. This also enables us to use 
a spectral expansion in the axial direction. So the only boundaries of concern are the 
transverse numerical boundaries. If we have good boundary conditions we can make the 
computational domain smaller and thus get a higher resolution and compute finer scales. 

Several investigators 9-14 have studied and used different types of non-reflecting 
boundary condtions. Engquist and Majda 12,13 study the wave equation and develop 
a perfectly non-reflecting boundary condition using pseudo-differential operators. But this 
equation is non-local in both space and time and is of little use for computational pur- 
poses as we would have to store data on the boundary from previous time-steps. But by 
approximating the operator they also develop a hierarchy of approximate NRBCs , with 
the increasing accuracy in passing obliquely incident waves. 

We use the first order non-reflecting boundary conditions for outflow at the transverse 
boundary, which can also be derived from physical considerations using 1-D Riemann 
Invariants. The outflow boundary condition is, 


& r ^ / Coo / 

P ~ PooCoo-U, - -—p 


at 


( 19 ) 


where R is radius of outer boundary, and the primed quantities represent the pertur- 
bations from mean flow and the oo quantities represent free stream or mean flow. The 
term on the right hand side is a correction for a cylindrical boundary. 

This boundary condition is supposed to ensure that there are no incoming waves from 
infinity. It works well for directly incident waves. We think that this will serve the purpose 
since we would have nearly cylindrical wavefronts incident on the boundary which is also 
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cylindrical. In addition we hope to mitigate the outgoing wave amplitudes by coarsening 
the mesh gradually as we approach the boundary as was done in the 2-D preliminary work. 

Now to implement this boundary condition in our scheme we have to transform the last 
set of equations, which represent the physics at the boundary, to characteristic variables 
at each substep of the Runge-Kutta algorithm. 

This is done using the following transformation matrix from Giles. 9 

/ c i\ / pL c Io 0 0 0 1\ 6a \ 

C2 Poo^oc'U'O Coo 0 0 0 StTIq 

c 3 = -ploCooUr 0 -Coo 0 1 Sm r (20) 

C4 Poo Coo 0 0 Coo 0 , StTI z 

V £5 ' ' PcXD^OO ^T 0 Cqo 0 1 / V Sp / 

where the C{ s are the characteristic variables and the 5-quantities are time incremental 
quantities. 

So we can modify the last five rows of the mass matrix and the right-hand-side using 
the above transformation. 

Equation 19 can be written as 

boo C'oo ^7' 5fX C oo (5 ! i l -|- Sp — “ —p 

2 Jr t 

So to implement the above b.c. we have to replace the third of the last five rows 
including the right-hand-side by the above equation since the third characteristic C3 has 
the right form for the b.c. 
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Test Cases 


A few simple test cases were run to validate the code and for diagnostic purposes. In 
all cases the desired steady state solutions were obtained to machine accuracy. All cases 
except the rotating flow case were run using quadratic B-splines. 

Case 1 : Uniform flow through a circular pipe 

Flow conditions are, 
uq 0, u r — 0, p — ~, a — T 
R = 1 


“> = l7<* 2 - r2 > 


Choose k = then, 


Uz = (1 - r 2 ) 

From the temperature equation (eqn 5) we have a balance, 

heat flux = - viscous dissipation 
or, 

^ V 2 T = ^ I- $ 


RePr Re 

(refer Appendix) 

Solving the above equation and applying the b.c’s 


dT, 

a7 lr=0 ~° 

T r =R — 1 


gives, 


T = ~^Pr( l-r 4 ) + l 
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This test case helps check the temperature equation and uniform flow in the z 


direction. 



Fig. 3 


Case 2 : Uniform flow through permeable boundaries 


Flow conditions are, u r = ucosO , uq = — usind , u z = 0, p = a = 1 
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Fig. 4 

This checks the flow jacobians and also non-axially symmetric flow. 

Case 3 : Uniform Rotating Flow 

Flow conditions are, 

ug = rQ, u r = 0, u z = 0, T = 1 
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The momentum equation reduces to, 


dp = puj 

dr r 
-yap = 1 

if p at outer radius is p 0 , 


■yr 2 fi 2 

p = p Q e 2 



Fig. 5 

This flow tests the azimuthal dependence. Here we had to switch to higher order B- 
splines since the variables are exponential functions and are better represented by higher 
order B-splines. 
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Case 4 : Spherically Radiating Flow 



Fig. 6 

This is to test the non-reflecting b.c’s. 

We consider a spherically radiating point source flow placed at the center of a cylin- 
drical domain. The solution is singular at the origin and is therefore picked up after a 
finite time and is treated as an initial value problem. This test case is very robust, since 
there exists an exact time-dependent solution with which we can compare our results. This 
test case not only tests the non-reflecting boundary conditions but also helps validate the 
entire code. 

The radiating bubble is considered as a perturbation to the mean flow which is, 

P = cr = 1, Ur = u z = U0 = 0 

The perturbation is taken as a point source solution of the linear wave equation which 
is written as a velocity potential $ = — Where x = \/r 2 + z 2 , c is sound speed, and 
t is the time. / can be any function. We have used a cubic function for /, 
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J^('-f) 3 (*i-(+f ) 3 c(t - tl ) < X < ct, 


<f> = 



The function is plotted below. 


otherwise 



Fig. 7 


The fundamental perturbations are expressed in terms of the velocity potential as, 


P 


- Poo 


dt 


o' — p' 




Ur dr 


Uc 


0 


Uz dz 


Plotted below are the density profiles as time progresses. The dotted profiles which 
almost overlap the solid lines are the exact solutions. Figure 8 plots density versus r at 
z = 0, while figure 9 shows density versus z at r — 0. Figure 10 shows a cross section 
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of the spherical wave in a r — z plane. From the plots we can see that the non-reflecting 
conditions work fairly well. The waves which are incident normally pass through the 
artificial boundaries with little or no reflection. In figures 8 and 10 the plots with t = t 3 
show a weak reflected residual wave. This is as expected since the non-reflecting conditions 
are exact for normally incident waves, but oblique waves would cause some reflection. These 
reflected wave amplitudes can be further reduced by stretching the grid near the boundary 
(not done here). 




Fig. 8 
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Fig. 10 

A more general and similar test case would be to shift the origin of the point source 
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off the z-axis. So the flow is not axis-symmetric anymore. 


Conclusions 


A Numerical code based on the spectral method of Moser 7 was developed in cylindrical 
coordinates using Fourier expansions in the azimuthal and streamwise direction and a 1-D 
B-spline basis representation in the radial direction. Time integration is carried out by a 
3 rd order Runge-Kutta time marching scheme. 

The code was validated against some test cases. Non-reflecting boundary conditions, 
used at the outer boundary, have also been satisfactorily tested. 
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Appendix 


The compressible Navier-Stokes equations in cylindrical coordinates 


dp 1 djrpv y) 1 djpug) d{pu z ) 

dt r dr r dd dz 


= 0 


(A.l) 


Du r 


V Dt 
Du e 


where, 


Dt 


+ 


_ 

r 

u r ue 


Du z 

) ~Dt 


dP_ dTrr 1 dTrd dr^ T rr - Tge 

dr dr r dd dz r 

\dP_ drer 1 dreg dr 6z r r9 

r d 9 dr r dO dz r 


dP dr zr 
dz dr 


1 dr z Q dr. 


+ r dO 


+ 


dz 


-A 

T 


(A.2) 


r = p 



l&iz. , 

du z | 9u r 
9r ' 9z 

1 T5f + ’•#■(“) 

2 (?^ + ^ - s v • fl ) 

9ug . 1 du z 

dz ^ r dd 

du z . du r 
dr ' dz 

9l/0 1 9u 2 

az r do 

2(^- — I V • 

Z ^ dz 3 V 




p is the viscosity, 


D d d ue d d 
Dt ~ dt +Ur dr + Vd 0 +Uz d^ 
_ _ 1 dru r 1 due du z 

r dr + r dO + dz 


dP 

dt 


+ V • Pu = —(7 — 1)PV • u + (7 — 1)V • q + (7 — 1)$ 


(A3) 


q is the heat flux, T is the temperature, k is heat conductivity, 


q = —kVT, T = 7 cr P 


and $ is viscous dissipation given by, 




1 due 


du. 


, 1 du z dug \ 2 ,du 

As? + 37) 


T + ) 2 + (i a * 

/ v nr' 


dr 
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r dug 


dr 
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